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O . 

, Non-Abelian plasma instabilities play a crucial role in the nonequilibrium dynamics of 

04 ' a weakly coupled quark-gluon plasma and they importantly modify the standard pertur- 

bative bottom-up thermalization scenario in heavy-ion collisions. Using the auxiliary-field 
formulation of the hard-loop effective theory, we study numerically the real time evolution 
of instabilities in an anisotropic collisionless Yang-Mills plasma expanding longitudinally 
in free streaming. In this first real-time lattice simulation we consider the most unstable 
modes, long-wavelength coherent color fields that are constant in transverse directions and 
which therefore are effectively l-|-l-dimensional in spacetime, except for the auxiliary fields 
which also depend on discretized momentum rapidity and transverse velocity components. 
We reproduce the semi-analytical results obtained previously for the Abelian regime and 
Oh| we determine the nonlinear effects which occur when the instabilities have grown such that 

^ ■ non-Abelian interactions become important. 

m ■ I. INTRODUCTION 

The experimental results obtained at the Relativistic Heavy Ion Collider RHIC [jj] and their 
I good agreement with hydrodynamical simulations with extremely early thermalization and a low 

shear viscosity jl, 0, 0] close to a conjectured lower quantum theoretical bound [5] are now widely 
. interpreted as evidence that the hypothetical quark-gluon matter produced at RHIC is strongly 

00 ! interacting and very far from a perturbatively accessible regime. Indeed, perturbative approaches 

■ like that of the original bottom-up thermalization scenario [a, 0, 0] do not seem to be able to come 

close to explaining the fast apparent thermalization. However, as pointed out first by Ref. 
the original bottom-up scenario is qualitatively changed by the inevitable presence of non-Abelian 



X 



' (chromo-Weibel) plasma instabilities [10|, llll, llj] in a weakly coupled quark-gluon plasma with 

^ - momentum-space anisotropy, although it is still an open theoretical question how the bottom-up 

scenario will have to be modified, even at asymptotically weak coupling and in the first stage of 
the bottom-up scenario 0, Q, m 0, Non-Abelian plasma instabilities have moreover been 



argued to importantly modify weak-coupling results on the shear viscosity to anomalously low val- 
ues (isl . Even if the quark-gluon matter produced at RHIC may be too close to the deconfinement 
phase transition for any extrapolations of weak-coupling results, it is clearly necessary to better 
understand the latter and how they differ from other approaches. Finally, it may be the case 
that the higher energies to be reached at upcoming heavy-ion collider experiments at the Large 
Hadron Collider (LHC) open the window to the specific collective phenomena of a weakly coupled 
quark-gluon plasma, such as non-Abelian plasma instabilities. 

In this paper we shall discuss only the theoretically clean situation at asymptotically weak 
coupling and the dynamical evolution of non-Abelian plasma instabilities in a collisionless plasma 
with long-wavelength color fields. Any amount of momentum anisotropy in the distribution of 
the (high-momentum) plasma particles leads to chromomagnetic instabilities, which in the weak- 



field situation are straightforward generalizations of the Abelian Weibel instabilities 19|] and whose 
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dis per sion laws have been worked out for specific cases of a stationary anisotropic plasma in Ref. [2C 



2ll. l22l. l23l|. In an Abelian plasma, the Weibel instabilities grow exponentially until they are large 
enough to modify the distribution of the hard particles and give rise to their fast isotropization. In a 
weakly coupled non- Abelian plasma, the situation is more complicated because the long-wavelength 
color fields have nonlinear self-interactions before they reach the size where fast isotropization 
occurs. The first numerical simulations 2j] of non- Abelian plasma instabilities using the systematic 
framework of the hard-loop effective theory 25, 2^, 27, 2^ have concentrated on the most unstable 
modes which are constant in the directions transverse to the direction of momentum anisotropy. 
It was found that such configurations experience a certain amount of Abelianization over domains 
of finite size when they enter the nonlinear regime, which allows them to continue an exponential 
growth out of the hard-loop regime, confirming essentially the conjecture of Ref. [29l] formed from 
numerical studies of a toy model which showed virtually complete Abelianization. In spacetime, 
the corresponding evolution equations are 1-1-1 dimensional, which in the hard-loop effective theory 
are coupled to auxiliary fields that depend on the three-dimensional velocity of the hard particles, 
so that in conventional plasma physics these simulations would be termed 1D-I-3V. Fully 3-1-1 
dimensional simulations (3D-I-3V) later showed however that more generic field configurations in a 
plasma with fixed (moderate) momentum space anisotropy do not continue to grow exponentially 



in the strong-field regime, but enter a linear-growth phase [30, l3l|] by the formation of a cascade 
which pumps the growing energy in the infrared modes into higher-momentum modes 15|, |32| ]. 
The recent simulations of Ref. [33] however found a continued exponential growth of initially small 
perturbations in the case of very strong momentum anisotropy. A very strong anisotropy (if not the 
requirement of initially small fluctuations ) is of particular interest for heavy- ion collision where 
in a weak coupling situation the longitudinal expansion makes longitudinal momenta of quarks and 
gluons much smaller than their transverse momenta. 

Recently, in Ref. [sj] the hard-loop effective theory for stationary anisotropic plasmas was 
extended to the case of a boost-invariant longitudinally expanding distribution of plasma particles, 
the hard-expanding-loop (HEL) effective theory. The essentially Abelian weak-fleld regime was 
worked out semi-analytically with the result that the counterplay of increasing anisotropy and 
decreasing plasma density lets Weibel instabilities grow exponentially in the square root of proper 
time, with more and more modes becoming unstable as time goes on, but each one experiencing 
a certain delay before growth kicks in. A similar behavior was previously found in numerical 
studies of initially small rapidity fluctuations in the so-called color glass condensate framework 
35I . 3^ . By matching the mass scales involved with the parameters of the saturation scenario 37 1 
the conclusion was drawn that LHC energies will be needed to allow for conditions where strong 
quark-gluon-plasma instabilities can develop from small initial rapidity fluctuations, leaving open 
however the issue of strong initial gauge flelds. 

In the present paper we begin the study of the evolution of genuinely non-Abelian plasma 
instabilities in a longitudinally expanding plasma by a lattice discretization of the HEL theory and 
1D-I-3V simulations. The latter have been found to give an upper limit of the full 3-1-1 dimensional 
evolution of more generic field configurations. The results of [33(] for strong anisotropy suggest 
that this upper limit may well be reached by 3-1-1 dimensional plasma instabilities that start out 
as small rapidity fluctuations (though not for those that are initially non-perturbatively large). 
3D-I-3V (as well as 2D-I-3V [38]) real-time lattice simulations of the HEL theory, which will be 
needed to address also initially strong rapidity fluctuations, will be the subject of follow-up work. 
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II. HARD-LOOP EFFECTIVE FIELD EQUATIONS FOR AN ANISOTROPICALLY 

EXPANDING NON-ABELIAN PLASMA 

For an ultrarelativistic plasma, a sufficiently small (gauge) coupling g introduces a hierarchy 
of scales, separating the hard momenta |p| = of plasma constituents from the "soft" scale 
~ g\/y\p\, where / is the typical hard particle occupation number (which may be different from 
order one in strongly nonequilibrium situations). The soft scale is associated with various screen- 
ing phenomena and the various branches of plasmon propagation. Ultrasoft scales ~ g'^f |p| are 
responsible for the damping of quasiparticles and, in or close to thermal equilibrium, for the non- 
perturbative screening of chromomagnetostatic fields. 

In an anisotropic plasma, the perturbatively accessible soft scale is also responsible for plasma 
instabilities, which constitute the dominant nonequilibrium effects at weak coupling: the associated 
rates are parametrically larger than any of the scattering processes, even though the latter are 
enhanced in a non-Abelian plasma.^ As long as the amplitude of the gauge fields A <C vTIpI) the 
evolution of the plasma instabilities is essentially Abelian and can be studied by a perturbative 
linear response analysis. For a stationary anisotropic plasma, the evolution is simply exponential 
in time. When the amplitude becomes nonperturbatively large, A > vTIpI; non-Abelian self- 
interactions of the gauge fields become important to leading order and require numerical evaluation, 
which as long as A <C \p\/g can be carried out consistently within the hard-loop effective field 
theory framework.^ In the latter, the hard particles are integrated out to produce a nonlocal and 
highly nonlinear effective action which can be written in terms of a compact integral representation 



43l . I44l . I45l | . This was initially obtained for the case of thermal equilibrium and has a straightforward 



generalization to the case of stationary momentum space anisotropy [26l . 1281 ]. It is of particular 



importance to numerical lattice studies that the corresponding effective field equations can be made 



local at the expense of introducing a continuous set of auxiliary fields 4611 w hich arise naturally 
when solving gauge covariant Boltzmann-Vlasov equations [isl . 47 . 48. l49l|. In the hard loop 
approximation, these auxiliary fields depend on the velocity vector of the hard particles whose 
hard momentum scale is integrated out. 



In Ref. 3j] this approach was extended to the case of a nonstationary plasma with a free 
streaming expanding distribution of hard particles, which we now review, filling in some details 
left out in Ref. [3^], before proceeding with numerical real-time lattice calculations. The latter 
allows us to follow the time evolution of plasma instabilities with initially small fields into the 
regime where non-Abelian self-interactions become important. The key difference to previous 
hard- loop simulations of non-Abelian plasma instabilities 24, 3^, 31, 33] is in the time-dependence 
of the (soft-scale) parameters which determine the growth rate of a given unstable mode and also 
which modes are unstable. 



A. Gauge-covariant Boltzmann-Vlasov equations in a nonstationary plasma 

Assuming a color neutral background distribution function /o(p,x, t) which satisfies 

vdfo{p,x,t) = 0, vf'=p^/p', (1) 



^ As we shall see below, for strongly anisotropic plasmas the relevant soft-scale parameters depend also importantly 
on the anisotropy parameter(s) hidden in /. 

^ For numerical simulations which take into account the backreaction of the soft fields on the hard particles that 
come into the play when A ~ \p\/g using a Boltzmann-Vlasov treatment see Refs. [sgI . |40| . IZH ]: for numerical 
simulations which include backreaction using a statistical classical field theory treatment see Refs. [ssl, [36.. .42.] . 
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the gauge covariant Boltzmann-Vlasov equations for colored perturbations 5 fa of an approximately 
collisionless plasma have the form 

V ■DSfa{p,^,t) = gv^Frdir^fo{p,^,t), (2) 

which have to be solved self-consistently with the non-Abelian Maxwell equations 

D,Fr = fa =9tRj 0^^^faiP^^^i)- (3) 

Here tR is a suitably normalized group factor, while the total number of degrees of freedom of the 
hard particles is taken care of by the normalization of the distribution function /q. 

In a stationary (but possibly anisotropic) plasma /o only depends on momenta, and ([T]) is sat- 
isfied trivially. Here we shall consider the generalization to a plasma which expands longitudinally, 
which should be a good approximation for the initial stage of a parton gas produced in a heavy ion 
collision as long as the transverse dimension of the system is sufficiently large. Assuming further- 
more boost invariance in rapidity fsB] and isotropy in the transverse directions, the unperturbed 
distribution function /o, being a Lorentz scalar, has the form 51. [s^] 

/o(p,x) = fo{p±,p\z,t) = fo{PL,p",T) (4) 

where the transformed longitudinal momentum is 



P 



''-i{p'-Pp°), P = z/t, ^ = t/T, T = V^^, (5) 



with = \ p\ + {p^y for ultrarelativistic (massless) particles. 



B. Comoving coordinates 

It is convenient to switch to comoving coordinates 

t = T cosh r], 13 = tanh 

z = rsinhr/, 7 = coshry, (6) 

i.e. a coordinate system with metric ds^ = dr^ — dx^ — r'^drf. We introduce the notation = 
= {t,x^ ,ri) with indices from the beginning of the Greek alphabet for these new 
coordinates. Note that in the latter the indices z, j, . . . are restricted to the two transverse spatial 
coordinates. 

In what follows we shall not deal with space-time covariant derivatives and Christoffel symbols, 
but write everything in terms of explicit derivatives. In particular the gauge covariant derivative 
always means^ Da = da — ig{Aa-, •]. Being a two form (where indices are naturally down), the field 
strength retains its usual form: Fap = daAp — dpAa — ig[Aa-, Ap\. The (non-Abelian) Maxwell 
equations do involve additional terms, but they can be written compactly as 



1 ^ . ^al3\ 1 



.Da{TF''^)^-Da Tg''^{T)g^'{T)F^& = f. (7) 

r r L J 

In addition to space-time rapidity rj, we also introduce momentum space rapidity y for the 
massless particles according to 

pf^ = p_L (cosh y, cos (/), sin (^,sinhy). (8) 



^ Recall that A'' = {(j),!) with 4-index up. Thus Ac = {Ar,-A'^, -A^, A^). 
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In comoving (tilde) coordinates, we then have 

= cosYirjp^ — smhrjp^ = p±cos\i{y — rf), (9) 
p^ = — Prj/r'^ = {cosh rjp^ — sinh'qp^)/ T 

= p'^/r = p_L sinh(y — ry)/r. (10) 

Instead of the Ught-Uke vector f ^ = p^^/p^ containing a unit 3- vector that was used in Eqs. ([1]) 
and ([2]), we shall define the new quantity 



p 

= — = ( cosh(y — ry), cos cp, sin (j), — sinh(y — r]) 
p_L \ T 

which is normalized so that it has a unit 2-vector in the transverse plane. 



1 



(11) 



C. Longitudinally expanding free streaming background solution 
Eq. ([1]) , involving space-time derivatives at fixed and p^ , can be rewritten as 



ip-d)fo 



0. 



y,PA 



Because 



p'^drPrjix) 



y,p± 



-p\ sinh(y — rf) cosh(7/ — rf) 

-p'^drjPrjix) 



this can be solved by /o(p,x, t) = foip±,p^ix)) = /o(p±, -p'^(a;)r(x)). 
In the following we shall use^ 



(12) 



(13) 



/o(p,x) = /iso 



pI + 



,p'^T , 



(14) 



which corresponds to local isotropy on the hypersurface r = riso, and increasingly oblate momentum 
space anisotropy at r > Tiso (but prolate anisotropy for r < Tiso)- Since a plasma description does 
not make sense at arbitrarily small times and so time evolution will have to start at a nonzero 
proper time the time Tiso may be entirely fictitious in the sense of pertaining to the pre-plasma 
(glasma 3^, HH]) phase. This will in fact be the case in the numerical simulations below, where we 
shall start already with oblate anisotropy by choosing Tiso < tq. 

In a comoving frame, the energy density and pressure components of the hard particle back- 
ground can be determined by evaluating T^^rt. = (2^)"'^ / d'^Px^y p°'p^ fo, which yields 



part. \ 



ppart. 



ppart. 



(r) 
(r) 



rpTT 

part. 



2 part. 



_rpri 

part.?; 



1 



arcsm 



T 



4(f2 - 1) 

1 



1 + 



2-1 

f2-2 



£: 



1 



; arcsm 



£: 



2(r2 



1) 



1 arcsin 



(15) 
(16) 
(17) 



Notice that it would be straightforward to relax the assumption of momentum-space isotropy in the transverse 
directions. 
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where Siso = '?part.(''"iso)) T = t/tiso and we have assumed r > 1. For r 1 we have 



ppart. 



The energy density follows from 



vr 



ppart. 



-8- f-3 



(18) 



part. 



Q ppart. . ppart. 



The particle distribution function (fT4l) has the same form as the one used in Refs. [20|, l2l|, l2J, |3l|] , 
but the anisotropy parameter ^ therein^ is now space-time dependent according to 



= (r/Tiso)' - 1, 



(19) 



and the normalization factor N{£^) of Ref. 2l|, |2J, |3l| is unity. 

The behavior ^ ~ at large r is a consequence of having a free-streaming background distri- 
bution. In a more realistic collisional plasma, ^ will have to grow slower than this. In the first 
stage of the original bottom-up scenario fis'], ignoring plasma instabilities, one would have had 
^ ~ r^/^. In Ref. 14 1 it was argued that plasma instabilities reduce the exponent to ^ ~ t'^^'^i 
whereas Ref. 17 1 recently presented arguments in favor of ^ ~ r^/'*. All these scenarios have ^ ^ 1, 
so below we shall concentrate on the case rjso < tq and thus high anisotropy for all r > tq, but in 
the idealized case of a collisionless free-streaming expansion. 



D. HEL effective field equations 

Transforming the gauge-covariant Vlasov equation to comoving coordinates one can write 

V -DSfX^ = 5y"F„"^af^)/o(p±,P,), (20) 

where the derivative on the left-hand side has to be taken at fixed as opposed to fixed p". On the 
right-hand side the derivative with respect to momenta is at fixed x, but the transformation from x 
to X does not depend on momenta anyway. However, in the following it will be important to write 
the right-hand side in terms of 9^^/o(p_L,p^) with index up so that this factor depends only on 

and and not additionally on r. This means in particular that p-d (5|^)/o)|p = p-d (c^^>)/o)|p = 0. 
Eq. ()20p can then be solved in terms of an auxiliary field W/3(x; 0, y) which satisfies 

V-DWp\^,^ = V^Fp^, (21) 

and 

5f{x-p) = -gWp{x;4>,y)d^^^Up±,Pr,). (22) 

The field VF/3(x; (j), y) is indeed analogous to the auxiliary field Wy{x; v) of the (static) hard-loop 
formalism [2^ because for a given space-time point it only depends on the 3-velocity of the hard 
particles, v = (cos </), sin 0, sinh y)/ coshy, and not on their energy p^ . Notice that only with index 
down its equation of motion (j2ip is formally the same as in the static situation. 



^ The anisotropy parameter 6 used in Ref. [ij is related to ^ by ^ ~ ^ 
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Expressed in terms of the auxiliary field W, the induced current in comoving coordinates reads 



(23) 



where for each ((/>, y) (i.e., fixed v) the scale pj_ (related to energy by p'^ = p^coshy) can be 
integrated out. 

With a distribution function that is even in p_L and pr^ as in ()14p , covariant current conservation 
can be verified without having to integrate partially with respect to p. (This proves to be helpful 
for the lattice discretization below, where all integrals will be replaced by discrete sums.) The 
current in ordinary coordinates is given by Eq. (j23p by dropping the tilde on and only. 
Starting from the first line of (f23l) . we can then use D-p = p- D = p-D, and (p ■ d)dfo/dpp\p = 
and finally (|2ip (with V replaced by p"). Changing the integration variables to p like in the second 
line of (f23]l we obtain 

. . fd'p^dp. 



(2^)' 2rJpi+p2/r2 



9p* dprj '^^'^ 

This vanishes already by symmetry when dfo/dp^ and dfo/dprj are odd functions in p* and p^, 
respectively. 

Specializing to the background distribution function (jl4p we have 



0, — COS 0, — sin (f), — 5- sinh(j/ — r/) j 



and we get 



1 + sinh^(?/ — rj) 



r = -rnl i j^^ ^ 1" (l + smh2(y " ^7)) >V(x; 0, y) , (25) 



(24) 



where 



= (cos </>, sin (j)) , Vrj = —T sinh(7/ — r\) , (26) 

and 

ml = -ghnl^^fUp)- (27) 



The mass parameter rriu equals the Debye mass at the (possibly fictitious because pre-plasma) 
time Tiso. 
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Because V ■ D commutes with the coefficients of Wi and appearing in the definition of W 
(in particular [V ■ D,Vr^] = 0, cf. p^ ). we do not need to evolve the components Wp separately, 
but only the combination W, which is governed by 

V-DW= (v'F.r + ^V^'Fr^r^ + V'V'^F,^ (28) 

V ^iso / V ^iso / 

For T — Tiso only Fq^-y (the electric field components in the comoving frame^) appear on the r.h.s., 
whereas for r 7^ Tiso magnetic fields come into the play, opening the door for magnetic instabilities. 

This single equation for W together with the Yang-Mills equations and the algebraic relation 
between j and W closes our equations of motion. To solve them, we adopt the comoving temporal 
gauge = and introduce canonical conjugate field momenta for the remaining gauge fields 
according to 

= TdrAi = -TdrA' = -Hi, (29) 

and 

= ^drA^ . (30) 

Notice that transverse (comoving) electric field components differ from 11* by a factor of r: 

= r-itf . (31) 

In contrast to most of the literature on the color glass condensate framework, we shall reserve the 
symbol E for the electric field and denote the canonical conjugate field momenta by 11. 
In terms of fields and conjugate momenta, the Yang- Mills equations take the form 

rd^W = jr, - DiF\ , (32) 
T-^drlii = f - DjF^' - Dr^F^' . (33) 



E. 1D-K3V equations 

A linear response analysis (appropriate for small gauge field amplitudes) shows that the most 
unstable modes of an anisotropic plasma are those whose wave vector is oriented along the direction 
of anisotropy. 

We therefore begin by considering only initial conditions and thus solutions which are constant 
in the transverse directions (i.e., neglecting transverse dynamics), SjA" = 0. Hence, -D* = —ig[A^, ■] 
and the Yang-Mills equations reduce to that of a 1-1-1 dimensional theory with A^ acting as adjoint 
scalar s. 

We then have 

-drUi = f + g^[A\i[A\A']] + ^D^^A' , (34) 
rdrll^ = jr^ + ig[A\D^A'], (35) 
as dynamical Yang- Mills equations, and 

Tf = D^W - ig[A\Ui] , (36) 



^ From here on we shall drop the tilde on the quantities in the comoving frame, which will be used exclusively in 
what follows. 
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as Gauss law constraint. 

The current j"' is a linear functional of W, given by Eq. (I25p as before, but the equation of 
motion for W, Eq. (p8]) . reduces to 



drWir,,;^,y) = ^^^^J^I?, ((l - .M^ - W 



r^MA^, W] + '-vm, - -'^^-^(y-^) u^ . (37) 

cosh(2/-??) r r4 

All fields here depend on the two remaining space-time variables r, rj, and the auxiliary adjoint- 
scalar field W additionally depends on the momentum space variables (f), y which parametrize the 
3- velocity in the colored fluctuations Sf"', cf. Eq. ([22]) . 

In the present paper we shall restrict our attention to this dimensionally reduced situation, 
which in conventional plasma literature would be referred to as 1D-I-3V, postponing the study of 
the more general 2D-I-3V and 3D-I-3V cases to future publications. 



III. LATTICE DISCRETIZATION AND NUMERICAL RESULTS 

A. Methods 



For a numerical evaluation of Eqs. (I34p - (l37p together with Eq. (I25p we discretize proper time 
starting with finite tq > and time step e. The space-time rapidity coordinate rj is made pe- 
riodic and discrete with Nfj points and (dimensionless) spacing a covering a rapidity interval 
(— A^^a/2, A^^a/2). The (matrix-valued) fields A^, , and W^^y are defined on the sites of the 
1-dimensional rapidity lattice, while the conjugate momenta n^,, Hy, and 11'' are defined on the 
temporal links. The gauge field Afj is replaced by the spatial link variable U = expigaA^i. 

The integration over the momentum-space variables (j) and V in Eq. (j25p has to be discretized 
such that covariant current conservation is preserved manifestly. When expressed in terms of (p 
and y integrals, the integrand in (j24p is either odd in y — or multiplied by sin0 or cos0. In 
order that discretization of y and (p respect manifest covariant current conservation, we thus need 
to respect reflection invariance in and y — The angular variable is made discrete with uniform 
spacing 2'ir/N^, but for y = y — rj we shall consider two possibilities. In method A we shall 
discretize the interval — A^^ ^ y ^ -^y uniformly with spacing 2Ay/(Ny + 1), and in method B we 
make the substitution y = atanh x and discretize the range — 1 -|- Ax < x < 1 — Ax with uniform 
spacing Ax = l/N^- Because of the rj dependence of the shifted variable y, the lattice equation of 
motion for the auxiliary fields W(f,^y that live on the y boundary have to be completed by boundary 
conditions for W in the y variable. For the W fields we do not impose periodicity, but instead take 
the Neumann condition dW/dy = at the y boundary. 

In Fig. [T] we show the evolution of the conjugate momentum Ila; oc cos{i'rj) in the case that the 
gauge group is taken to be Abelian U(l). The system is initialized with a single Abelian U(l) mode 
with only Hx initialized with rapidity wave number u = 167r/5 = 10.053. . . in order to facilitate 



comparisons with semi-analytic results obtained in an earlier work [3j] where, for the Abelian case, 
the equations of motion for the W field have been solved in terms of integro-differential equations. 
Fig. [1] compares a semi-analytic result obtained from the latter with results obtained using the two 
different methods of discretization described above and detailed in Apps.[A]and[Bl As can be seen 
from this figure both numerical discretizations reliably reproduce the Abelian U(l) semi-analytic 
result. In the inset we compare the relative error defined as the difference of the time evolution 
obtained from methods A or B with the semi-analytic result over the sum (relative percentage 
error). As can be seen from this inset method B seems to perform better at late times so unless 
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FIG. 1: Proper-time evolution of the canonical field momentum Tlxir] = 0) of a single Abelian mode with 
rapidity wave number i/ = 10.053. Solid line is the semi-analytic result of Ref. [34| and the dot-dashed 
lines are the results obtained from our 1+1 numerical solutions using two different methods (A and B) for 
discretizing the shifted momentum-space rapidity y ~ y ~ -q. Inset shows relative error of the two methods. 
Run was made using n^o = 0.1, tq = 1.0, ruD = 10, o = 0.0025, e = 0.001, iV^ = 250, = 8, and (Method 
A) Ny = 1000, (Method B) = 1000. 



otherwise indicated all final results presented will be using method B. However, in practice, we 
have made runs comparing the predictions of methods A and B in all cases and find that there is 
very little difference between the results obtained with the two methods. 



B. Single-mode results 

In Fig. [2] we show results of a simulation of a single SU(2) mode with rapidity wave number 
u = 10.053 (same mode as Fig. [T]but now also with the color direction rotating with period 27r/z/ 
in space-time rapidity rf). In Fig. [2^ we show the proper-time evolution of the magnetic, electric, 
and total field energy densities in units where tq = 1 and, following Ref. |3Q], scaled with a factor 
of T. Because the energy in the hard particles is dropping proportional to t~^, this corresponds to 
giving the various soft energy densities in terms of the hard energy density (times a parametrically 
small number ~ since the hard energy density is assumed to be much larger than the soft ones 
in order that the hard-loop approximation be applicable.) 

The various components of the (soft) field energy density are defined by 

£ = £t + £l = £bt + + + £et 



tr 



t-2f2 +r-2nf + +(n^)2 . (38) 



Because of the expansion of the system, the total energy density £ is not conserved, even when the 



induced current (|25p is identi cally zero. In this case the time dependence is governed by the fact 



that the Hamiltonian density [36] TL = t£ satisfies 



and therefore 



|:« = ^W = £.-f., (39) 



^^\m = -'^^r\jm- (40) 
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FIG. 2: Results from a run with a single non-Abelian mode with v — 10.053. In the top panel (a) we 
show the proper-time dependence of the chromo-field energy densities. In the lower panel (b) we show the 
longitudinal and transverse pressures along with our numerical Gauss law violation. Run was made using 
Tiso = 0.1, To = 1.0, mo = 10, a = 0.0025, e = 0.001, = 500, = 100, and = 100. 



In the presence of a plasma of hard particles and thus nonvanishing induced current j we define 
the net energy gain rate by 

_ d£ 2 

Energy Gain = + , (41) 

which in the plots showing the energy densities is included as the dotted line marked "Gain Rate" . 
The latter gives the rate of energy transfer from the free-streaming hard particles into the collective 
chromo- fields. As can be seen from Fig. [2] for SU(2) the single mode evolution is quite compli- 
cated with all field components being dynamically generated; however, at late times transverse 
chromoelectric and chromomagnetic fields exponentially dominate. 

In Fig. [2)3 we plot the longitudinal and transverse field pressures generated during the system's 
dynamical evolution. These are obtained from [s^ 



Pl = £t-£l, 
Pt = £l, 



(42) 
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FIG. 3: Results from a run with a single non-Abelian mode with rapidity wave number v = 10.053 in which 
we have decoupled the hard particle currents (j — 0) so that we are simply solving the Yang-Mills equations 
in the expanding metric. In the left panel (a) we show the proper-time dependence of the chromo-field 
energy densities. In the right panel (b) we show the longitudinal and transverse pressures. Gauss law is 
obeyed exactly by our algorithm in this case. Run was made using Tjso = 0.1, tq = 1.0, ttid = 10, a = 0.0025, 
e = 0.001, and = 500. 



where as before Et is the sum of the energy density coming from transverse electric and magnetic 
fields and Ei, is the sum of the energy density coming from longitudinal electric and magnetic 
fields. As shown in Fig. ^Bp the system generates both longitudinal and transverse pressures. At 
short times (t/tq ~ 5-6) for this single mode evolution we find that the longitudinal pressure 
becomes momentarily negative; however, at late times the effect of the chromo-field instability is 
to generate exponentially large longitudinal field pressure, whereas the longitudinal pressure of the 
(free-streaming) particles drops according to tP'j^^^' ~ r~^. 

Also shown in Fig. [2)3 is our measure of violation of Gauss law which is determined by evaluating 
the r-component of the equations of motion as detailed in Eqs. (I36p and (IA13j) . As can be seen from 
this figure although our violation of the Gauss law constraint grows with time, it is numerically 
under control and always orders of magnitude below the field energy density. The amount of 
violation can be systematically reduced by taking finer lattices in r/ and velocity space. We have 
found that our results for the time evolution of the energy densities, pressures, etc. remain the 
same as our numerical Gauss law violation is reduced giving us confidence in our algorithm. As a 
general rule we have always terminated our runs when the Gauss law violation becomes of order 
one. 

For comparison in Fig. [3] we show the evolution of the field energy densities in the case of pure 
Yang- Mills evolution. This is obtained by decoupling the free-streaming particle currents by setting 

to zero in the field equations of motion. From Fig. [3^ we see that in the case of pure Yang- Mills 
evolution the field energy density decreases over the entire time interval shown. The "Gain Rate" 
control variable is approximately zero and shows the level of discretization errors. In addition 
we see that although both longitudinal and transverse pressures are generated they are of much 
smaller magnitude than those generated when the free-streaming particle currents are coupled into 
the Yang-Mills equations. Therefore, we have demonstrated that coupling in the particle currents 
generates qualitatively different field dynamics. 
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FIG. 4: Results from non-Abelian run initialized with a random superposition of discrete electric modes 
(cutoff white noise). In the top panel (a) we show the proper-time dependence of the chromo-field energy 
densities and the energy gain rate (PT|) times an extra factor of tq. In the middle panel (b) we show the 
longitudinal and transverse pressures along with our numerical Gauss law violation. In bottom panel (c) we 
show the correlations ■fyib], and C[j]. Run was made using Tjso = 0.1, tq = 1.0, uid = 10, a — 0.03, 
= 20, a = 0.0025, e = 0.00025, TV,, = 1000, = 100, and = 100. 
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FIG. 5: Results from non-Abelian run initialized with FGM initial conditions. In the top panel (a) we show 
the proper-time dependence of the chromo-field energy densities and the energy gain rate (|4ip times an 
extra factor of tq. In the middle panel (b) we show the longitudinal and transverse pressures along with our 
numerical Gauss law violation. In bottom panel (c) we show the correlations Cbli ^^id C\]\. Run was 

made using riso = 0.1, tq = 1.0, = 10, cr = 0.05, = 20, a = 0.005, e = 0.0005, TV,, = 500, N.^ = 200, 
and = 200. 
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C. Initial non-Abelian noise and partial Abelianization 

In Fig. |4]we show results obtained from an SU(2) non-Abelian run in which the initial condition 
is taken to be a random superposition of discrete transverse electric field modes (both Ux and Uy) 
with an ultraviolet cutoff Ai, = 20 in space-time rapidity wave number i/. The amplitude for each 
mode is chosen from a Gaussian probability distribution centered at zero with standard deviation 
a = 0.03. As can be seen from Fig. [4^ the system very quickly generates chromomagnetic fields 
whereas during the early times energy is equally distributed between transverse chromoelectric and 
chromomagnetic fields. Longitudinal field energies which vanish initially grow exponentially with 
a rate about twice of that of the transverse fields, but almost saturate when the nonlinear regime 
is reached. During the initial growth phase as well as in the deep nonlinear regime, the energy 
density is exponentially dominated by transverse chromomagnetic fields. This again translates into 
the generation of exponentially large longitudinal pressure as shown in Fig. [Ha. 

In Fig. we plot various measures of the Abelianization and (color) correlations of the chromo- 



fields. Following Ref. 29|, |3l| we define a measure of the "Abelianness" of the field configurations 
through 

If the field configurations are Abelian (aligned in one color direction) then this quantity vanishes 
because of the commutator in the numerator. 

In order to further study the color correlations of the chromo-fields in spatial rapidity, r], we 
define 

(P\ = -1 dr, tr {(i[ji{r] + 0,U{ri + ^.r^)3jiv)]?] 
^^^^^ 2N, io tr{j|(r? + 0}tr{j2(r?)} ' ^ 

where U{rj',rj) is the adjoint-representation parallel transport from tj to rj' . When colors are com- 
pletely uncorrelated over a distance ^, this quantity equals unity; if they point in the same direction, 
this quantity vanishes. Following Ref. 2^, 31 1 we define the "Abelianization correlation length" 



as the smallest distance where XA is larger than 1/2, 

U[j]= min (e). (45) 

Xa(0>1/2 

This we compare with a general correlation length, which does not focus on color, defined 
through the gauge invariant function 

Jo" dv^T^{ji{v)ji{v)} 

This function now vanishes when fields are uncorrelated over a distance ^, and it is normalized 
such that x(0) = 1- We thus define the general correlation length through 

m= min (e). (47) 

x{0<l/2 

Fig- Hb shows that the system becomes Abelianized with large color correlation length, 
when the fields have grown such that nonlinear self-interactions become important. occasionally 
even shoots up to the size of the space-time rapidity lattice (2.5 in this case) before settling to 
oscillations around rapidities ~ 0.3. (The indication of some late-time growth of is presumably 
spurious, since it is accompanied with the onset of a rapid growth of the Gauss law violation control 
parameter.) Although we show the output of only one run here the behavior shown is generic for 
all random seeds we have studied. 
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D. Color-Glass-Condensate-inspired initial conditions 

In Fig.[5]we show results obtained by using initial seed fields which reflect the spectral properties 
obtained by Fukushima, Gelis, and McLerran (FGM) within the Color-Glass-Condensate (CGC) 



framework 5j]. We use again a random superposition of modes, but now involving already initially 



both chromoelectric and chromomagnetic transverse fields with a spectrum^ 

\Ui{iy)\r=ro = 

\A\v)\r=r^ = al^, (48) 

for all space-time rapidity wave numbers v < that are allowed by the periodic boundary 
conditions of our finite t] lattice, excluding however = 0. The phases of each color component 
of these modes is taken at random, and we have used a small value o = 0.05, corresponding to 



initially weak fields. In accordance with Ref. 5j], the longitudinal magnetic field is set to zero 



initially through A^|T-=ro = 0, but the non-Abelian Gauss law leads to nonvanishing longitudinal 
chromoelectric fields even though = initially. 

With = initially, the Gauss law constraint in the 1+1-dimensional setting gives 

d^W = ig[A\Il,]. (49) 

Having populated the transverse field modes according to Eq. (j48p . we solve the lattice version of 
Eq. (I49[) to determine the longitudinal electric field = 11* /r. 

However, in contrast to the simpler initial conditions used above, this presents a problem with 
the periodicity of our rj lattice, since the solution thus obtained does not share the periodicity of 
all other fields, leading to a Gauss law violation at the boundary in the form of a mismatch of 
H''. This initially small violation however quickly grows and cannot be tolerated. We have solved 
this problem by singling out the lowest lying mode of and to calculate its contribution to the 
mismatch of H''. By elementary linear algebra we determine how to rescale the color components 
of this one mode such that the mismatch is eliminated, but this rescaling is only accepted when 
the total amplitude of this mode does not get modified by more than 50%. If this is not the case, 
a different set of random numbers for the phases of all transverse color fields is generated and the 
procedure repeated until a configuration is found where the amplitude of the lowest lying mode of 
Ax is not too far from the starting point ()48p . 

The results shown in Fig. [5] are qualitatively similar to those obtained with a random superpo- 
sition of purely electric modes initial condition. As can be seen from Fig. [5ti at early times there is 
equal partitioning between chromoelectric and chromomagnetic fields which both initially decrease 
and then begin to grow exponentially with transverse chromomagnetic fields dominating for nearly 
the entire run. At t/tq ~ 13 there is a non-Abelian "bounce" when the longitudinal field compo- 
nents become on the same order of magnitude as the transverse ones; however, beyond this point 
in time the transverse field components again dominate. In Fig. [Sb we see that the field pressures 
which are generated are also similar to those obtained with a random discrete Fourier spectrum 
with the system generating an exponentially large longitudinal pressure due to the chromo-Weibel 
instability. In Fig. [5}: the behavior of the Abelianization measure, ^[j], and correlation lengths are 
again similar to the random discrete Fourier spectrum initial conditions showing an Abelianization 
of the fields and large color correlation length at late times. This demonstrates that the qualitative 



^ The spectrum of fluctuation derived in Ref. [sj] of course has also modes which are not constant in the transverse 
coordinates, but in our present framework we have to restrict ourselves to modes which are effectively 1+1- 
dimensional. 
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FIG. 6: Field energy density results from an Abelian run initialized with FGM initial conditions. Transverse 
fields and Gauss law violation are zero to within machine precision. The field pressure is purely longitudinal 
and coincides with the total field energy density. Run was made using Tiso — 0.1, tq = 1.0, = 3.585, 
a = 0.01732, = 20, a = 0.01, e = 0.001, iV^ = 500, = 100, and = 4. 



features of the time evolution of the instability induced fields are independent of the details of the 
initial condition.^ 

In Figs. [6] and [7] we show results obtained using FGM initial conditions with a smaller Debye 
mass corresponding to the estimates of the "gluon liberation factor" c obtained from the color-glass- 
condensate picture [sH, (see Appendix [C] for details). In Fig. [6] we show the results obtained 
from an Abelian run in which all fields were constrained to initially point in the same direction in 
color space and in Fig. [7|we show the results of a non- Abelian SU(2) run. As can be seen from both 
figures the primary effect of lowering m^) is to slow down the growth of the chromo- fields; however, 
besides this "stretching" of the time axis there is little qualitative difference between the larger ttid 
run (Fig. [5]) and this case. Fig. [71 We still observe domination by transverse chromo- fields, which 
now have larger color correlation length, and generate exponentially large longitudinal pressure. 

In Fig. [8]we compare six different non- Abelian SU(2) runs with FGM initial conditions in which 
we have taken different values for the spectral cutoff in rapidity wave number, A,y, imposed on 
the FGM initial condition. As can be seen from this figure for fixed initial energy density the 
effect of increasing Ky is to delay the onset of exponential growth of the chromo- fields. This is to 
be expected since for fixed energy density the occupation number of the lowest v modes must be 
decreased as is increased, and higher modes have a larger delay, as already found in the Abelian 
case studied in Ref. [s^]. In fact, the amplitude of the low- momentum modes must be decreased 
rapidly since the high-momentum modes dominate the energy density. In Fig. [8j) we show a fit to 
the "time to return", tj?, of the scaled energy density ttq£, i.e. the time it takes the instability 
to compensate for the initial decay of the soft fields caused by the system's expansion. Fitting 
this time (in units of tq) by a power-law a{A,y)'^ we find a = 13.46 it 0.01 and b = 0.26 it 0.01. 
The exponent b is consistent with being 1/4. The coefficient a depends on the Debye mass and 
decreases as ijid increases. 

In Fig. [9] we compare the rapidity (z^) spectrum obtained by Fourier transforming the trace of 
the conjugate field momenta, tr(n'^) = tr(n? -|- r^(n'')^), in order to gain more understanding of 
the momentum space dynamics of the fields in our simulations. In the left panel. Fig. [9ti, we show 



Of course, by this we mean any reasonable initial condition. Choosing, for example, an initial condition which 
only had very high frequency modes would greatly delay the onset of instability driven growth of the fields. 
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FIG. 7: Results from non-Abelian run initialized with FGM initial conditions. In the top panel (a) we show 
the proper-time dependence of the chromo-field energy densities and the energy gain rate (j4ip times an 
extra factor of tq. In the middle panel (b) we show the longitudinal and transverse pressures along with our 
numerical Gauss law violation. In bottom panel (c) we show the correlations ^a^], ^^id C[j]. Run was 
made using Ti^o = 0.1, tq = 1.0, mo = 3.585, a = 0.01, = 20, a = 0.01, e = 0.001, = 500, = 200, 
and = 100. 
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FIG. 8: Left panel (a) total field energy density results from a non-Abelian run initialized with FGM initial 
conditions and different UV cutoffs, Aj^ S {5, 7, 10, 20, 50, 100}, imposed on the initial spectrum. Right panel 
(b) shows the "time to return" t/j/tq, defined as the point at which the field energy density has returned 
to its initial value, as a function of the FGM spectral cutoff Ki, on a log-log plot. Blue line shows a fit to 
a power law tr = aTQ{K^Y' . Runs were made using Tiso = 0.1, tq = 1.0, tod — 3.585, a = 0.06, a = 0.005, 
e = 0.0025, Nn — 1000, A^j, = 800, and A^ — 50. For this figure discretization method A was used. 




mode number (v) mode number (y) 

FIG. 9: Fourier spectrum of the color-traced conjugate field momenta, tr(n^), obtained from (left) Abelian 
and (right) non-Abelian runs with FGM initial conditions. The lowest (bold green) line indicates the starting 
spectrum and the uppermost (bold red) line indicates the final spectrum. In the right panel the bold blue 
line indicates the "non-Abelian point" at t/tq ~ 55 when all field components become approximately the 
same order of magnitude. The Abelian and non-Abelian spectra were obtained by analyzing the currents 
produced during the runs shown in Figs. [6] and El respectively. 

the spectrum resulting from analysis of the induced current from the Abelian run shown in Fig. [H 
In the right panel, Fig. [9)3, we show the spectrum resulting from analysis of the induced current 
from the non-Abelian run shown in Fig. [71 The lowest (bold green) line indicates the starting 
spectrum, the bold blue line indicates the "non-Abelian point" at which all field components 
become approximately equal in magnitude, and the uppermost (bold red) line shows the final 
spectrum obtained in our simulations. As can be seen from this figure there is a stark qualitative 
difference between the Abelian and non-Abelian spectra with the former maintaining the spectral 
cutoff imposed on the initial condition and the latter "cascading" energy to higher and higher 
momentum modes starting already at very early times. This is similar to earlier results for the 
spectra induced by instability growth [15]. Surprisingly, in Fig. [9)3 one sees that at the "non-Abelian 
point" indicated by the bold blue line that the low frequency modes have generated a quasi-thermal 
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(Boltzmann) distribution up to z/ ~ 80. In fact, the development of the quasi-thermal distribution 
begins at very early times and one can associate a temperature with the system by fitting the 
low-i^ spectra with exponential fits from rather early times. Similar spectra are generated when 
one measures tr(A^) which also allows one to define a kind of magnetic temperature from that 
observable as well. 



IV. DISCUSSION, CONCLUSIONS, AND OUTLOOK 

In this paper we have performed the first numerical study of non-Abelian plasma instabilities 
in a nonstationary, longitudinally expanding system within the framework of discretized hard loop 
theory^, extending the semi-analytical results of [sJ] for the weak- field, Abelian regime. We have 
worked out the case of the most unstable modes which are constant modes in the transverse direc- 
tion, making the dynamics 1+1-dimensional in configuration space (while momentum space remains 
3-dimensional) . Starting with only small rapidity fluctuations, we found that the exponential (in 
-y/r [3?]) growth in the Abelian (weak- field) phase is only mildly weakened when nonlinearities 
through non-Abelian self-interactions of the collective fields set in, and this is associated with sig- 
nificant degree of Abelianization in finite domains in the nonlinear regime. This is quite similar 



to what was observed in the 1D+3V simulations in a stationary anisotropic plasma [2j] and it 
remains to be seen what full 3D+3V simulations will give. However, it is quite plausible that the 
1D+3V results already capture the behavior of the more generic 3D+3V simulations, because it 
was recently observed [33.] that for extreme anisotropies a saturation of the growth as was found 



in Refs. [30|, |3l|] at moderate anisotropies will occur only at correspondingly extreme values of the 
fields, if at all. Indeed, our simulations start out with strong anisotropy of the particle distribution, 
which rapidly grows with increasing time according to Eq. (jl9p . 

In our simulations we have found that in the non-Abelian case the growing unstable modes tend 
towards a quasi-thermal spectrum (Fig. [Db) and they produce mainly longitudinal field pressure, 
which grows exponentially, thereby realizing a bottom-up isotropization scenario in which the soft 
modes make up for the strongly decaying longitudinal particle pressure, which goes like The 
transverse particle pressure, which according to the CGC picture is approximately thermal by 
itself, is decaying like 1/r. In the hard-(expanding)-loop theory which we have considered, we can 
of course only trust the beginning of this scenario, since the backreaction of the collective fields 
on the hard particle background is neglected, and it is in fact the reservoir of energy in the hard 
particle background that is feeding the growth of the soft modes, which has to stop before the 
energy in the latter becomes comparable with the former. 

In Fig. [TOk we have reproduced the results of Fig. [7)3 for the field pressures obtained by choosing 
dimensionful parameters motivated by the CGC scenario as described in App.Oand compared also 
with the particle pressures that follow from this matching. Notice that all quantities are multiplied 
by T so that the decaying transverse particle pressure is represented by an approximately horizontal 
line. 

The time scale tq ~ Qj^ can be roughly identified with 1 — 1.5 and 3 GeV for RHIC and LHC 
experiments, respectively, where the plasma lifetimes are probably less than 5 fm/c ~ 25 — 35ro 



for RHIC, and probably much larger than 7 fm/c ~ IOOtq for the LHC [571]. Defining an effective 
growth rate of the longitudinal pressure by 

7L = ^ln(rr3PL,ficid), (50) 



Closely related instabilities have been found before numerically in the color-glass condensate framework in Ref. [35 
[3^ . where the role of plasma particles is played by high-momentum modes of the Yang-Mills field. 
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FIG. 10: Top panel (a) shows a comparison of particle and field pressures generated during a typical run 
which uses CGC-inspired "FGM" initial conditions. Simulation parameters are the same as shown in Fig. [7] 
and the field pressures are the same as in panel (b) of that figure. The particle pressures are obtained by 
evaluating the expressions given in Eqs. and p?]) . In order to fix the initial energy densities we use the 
scheme detailed in Appendix [Cl assuming as = 0.3 which gives £{to) ~ 0.3 Q^. In the lower panel (b) we 
show the growth rates for the field pressures in units of Tq"^. 



we find for the example provided in Fig. [10] a maximal value of about 7l ~ 0.2 Tq^ in the weak-field 
regime for r > 20 tq, and about 0.1 Tq"^ in the strong-field (non-Abelian) regime r > 70 tq. This 
corresponds to minimum characteristic time scales of 

. _i r0.7-lfm/c (RHIC) 
mm 7r ~ s , . „s (51) 

\ 0.3fm/c (LHC) ^ ^ 

in the weak-field regime, and twice that in the strong-field regime. This agrees roughly with the 
pre-isotropization values obtained in Ref. [42] from classical-statistical lattice gauge theory. ^'^ 

However, at least for the case of initially small rapidity fluctuations which we have considered 
here, there is a delay of the onset of plasma instabilities caused by the expansion which appears 



The higher growth rate of the transverse field pressure, which is due to non-Abelian self-interactions of the chromo- 
fields (it vanishes in the Abelian case), is what Ref. 42] would call a "secondary" instability. 
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uncomfortably large for RHIC energies, even if one chooses smaller spectral cutoffs in the initial 
fluctuations which somewhat reduce this delay (cf. Fig. [8]). 

Still, for the LHC our results suggest that plasma instabilities like those studied here will be 
an important phenomenon, in particular if LHC energies make contact to a more weakly coupled 
quark-gluon plasma as suggested for instance by the analysis of Ref. [s^. The comparison of 
particle and field pressure in Fig. [TOk indicates upper limits for an isotropization point, which are 
however strongly dependent on the initial strength of the rapidity fluctuations. Larger seed fields 
will correspondingly lower this point. However, experience from simulations of non-Abelian plasma 



instabilities in the stationary anisotropic case [30|, l31|, IS^l lets us expect that full 3D+3V studies 
(or at least 2D+3V ones [So]) are required to analyse truly strong initial fields. This will be the 
subject of follow-up work. 
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APPENDIX A: LATTICE DISCRETIZATION METHOD A 

The one-dimensional situation studied herein assumes that fields vary only in the ry-direction. 
We then have transverse adjoint scalar fields withi = 1, 2 and auxiliary fields W^^y which all are 
defined on the sites s of a periodic spatial rapidity lattice with (dimensionless) lattice spacing a. 
The conjugate momenta Hj live on the temporal links, while the conjugate momentum H'' of the 
gauge field A^^, which appears only in the form of a parallel transporter U^_^_i = exp (igaAfj^s), will 
be treated as located on the timelike plaquette on top of the link between site s and s + 1. 

Apart from U, all of these fields are represented by Nc x Nc traceless Hermitian matrices which 
for SU(2) reduce to the 2x2 Pauli matrices. Although we are going to make explicit all occurrences 
of the coupling g, in practice we have taken g = 1 through a rescaling of the fields. 

Covariant derivatives are defined in three versions: left- and right-covariant, 

A'^-U,_^A'^__,UI, Ul A'^^,U^^,-A'^ 
D^A': = i ^ , Df;A^ ^ , (Al) 

and symmetric, 

D'^^{D^ + D^)/2. (A2) 

The second-order is given by 

Z)2 ^ (D^ - D^)/a (A3) 

and is automatically symmetric. 

In method A, the auxiliary field W(r, r/; (f), y) of the continuum theory is modelled by a large 
number of fields Ws;0,g with y = y — rj discretized with Ny points in the interval (—Ay, Ay) and N^j, 
points for < (j) < 2tt. Additionally, we can absorb all or part of the denominator appearing in 
Eq. psp for the induced current by writing 

m;^,yi'^) = r\r,y)Ws;^,yiT) (A4) 
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with 



/(r,y)= l + ^sinh2(y) , 



(A5) 



with A a number between and 2. This does not produce extra terms in the equation of motion 
for W because 



5. + 



/ = o. 



(A6) 



The induced current (|25p is obtained from the auxihary fields (which at r = tq = 1 are taken 
to vanish) according to 



(A7) 



with defined in pT]). 

The equations of motion of the various fields are then solved numerically by a leapfrog procedure. 
The first step is to calculate the conjugate momenta from 



1 



n,,,(r + -) = n,,,(r - -) + e ( tjI + -D^Al + Tg'i[Ai, i[Ai,Al]] 



n^(r + -) = n^(r - -) + 6 i^--{fj + c/]^.i,\if/,+i ) + ^[a:, I^^^A 

The second step is to update the fields according to 

Al{T + e) = 4(r) + 6(r+l)-X.(r + |), 
f/,+ i(r + 6) = exp(i5ea(r+i)n^(r+^))c/^^i(r) 



^5 r »i rM At- 



and the auxiliary fields W according to 



y^s;4>,y{r + e) = W.;<^,g(r - e) + 2e{/(r, yj-^C - 5 



cosh(y) 



tanh(y) 



with 



sm 



Hy) f. 

4r2„ I 



r + - 



n?(r + |) + f/._in^.i(r + |)c/]_. 
n^(^-^) + f^.-inLi(r-^)c/]_, 



The Gauss law constraint is checked by evaluating 



(A8) 

(A9) 
(AlO) 



} (All) 



(A12) 



1^ E ( ^D'^^'si^ +'-)- ^[A,(r), n,,,(r + '-)] - fir 



2' T 



(A13) 



and then taking a square root to obtain the results shown reported in the main body of the text. 
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APPENDIX B: LATTICE DISCRETIZATION METHOD B 

In method B, the shifted momentum space rapidity y = y — r] is not discretized with uniform 
spacing in y but in a velocity-hke variable x, — 1 < x < 1, defined by 



y = atanh(3;) , dy 



1 



zdx. 



(Bl) 



Compared to method A, this has the effect of giving more lattice points around y = 0, where the 
VV functions are typically sharply peaked. 

With sinh^(y) = — x^) and cosh^(y) = 1/(1 — x^) this leads to 



Ws;0,x(r + e) = Ws;4>,x{r -e) + 2e\ /{t, x)-' C - g{l - x^)^ 



+ 



/(r, x)'^ B - {Df^ - (1 - x^)d^)Ws;4>A^) 



(B2) 



with (choosing now A = 2) 



and 



/(r,x)^(l-x2) 2(l + (_-l 



X 



(B3) 



T^X 



n^(T-|) + c/._in^-i(r-|)c/]_ 



The currents are then given by 

f-- 



r-2TT r-1 y _ 

ml) J d(t> J dx{l-x^y^W, 
f = -m% del) dxv' {1- x'^y^W , 
d(l) dxx (1 - x^) 2 W , 

J-1 



f = -ml) 



where the integrations over x and 4> are replaced by uniformly spaced discrete sums. 
All other lattice equations of motion are as in method A. 



(B4) 



(B5) 



APPENDIX C: MATCHING TO CGC PARAMETERS 

For fixing the dimensionful parameters of our numerical simulation in a way that makes contact 
with heavy-ion physics, we proceed as in Ref. [sS] and refer to the Color-Glass-Condensate frame- 



work [371 . |59| and take as starting time for the plasma phase tq — , where Qs is the so-called 



saturation scale. 
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In order to determine the only other dimensionful parameter in our HEL effective field equations, 
the Debye mass mo aX the (fictitious because pre-plasma) time tjso, we assume a squashed Bose- 
Einstein distribution function for the hard particle distribution function (|14p through /iso(p) = 
M(2Ng)/{e^^'^ — 1) where Ng = — 1 is the number of gluons and M a normalization that is 
adjusted such that at r = tq the hard-gluon density of CGC estimates is matched. Since the 
expansion is by assumption purely longitudinal, T is a constant transverse temperature, and it has 
indeed been found in CGC calculations that the gluon distribution is ap pro ximately thermal in the 



transverse directions, with T = Qs/d and d ^ ~ 0.47 according to Ref. 37|. The normalization Af 



can then be fixed by following Ref. [60], who write the initial hard-gluon density as 

niro) = c (CI) 

where c is the gluon liberation factor, for which different estimates can be extracted from the 
literature. 



According to Ref. [6Cf], the numerical CGC simulations of Ref. [6ll . l62l | correspond to c ~ 0.5, 
while an approximate analytical calculation by Kovchegov [55i] gave c = 21n2 ~ 1.386. We adopted 
this higher value for the numerical simulations in Figs. [6]l9l which is the more optimistic one from 
the point of view of plasma instabilities and which is actually not far from the most recent numerical 



result c ~ 1.1 by Lappi [56|] 



With Tiso remaining a free parameter which determines how anisotropic the gluon distribution 
is at To, the normalization M is now fixed by 

n(ro)^ = n(ri,o) = ^AfNgT\ (C2) 

Tiso 

For a purely gluonic plasma, the isotropic Debye mass is given by 

mj:,[Tiso)=N , (C3) 



which together leads to 



ml{riso)rl{Qsro)-' = - 1.285^, (C4) 

dCCj) Tiso Tiso 



when c = 2 In 2 and Nc = 3. We adopt this value for our simulations where = 2, since in previous 
studies of the stationary anisotropic situation little difference was found between the SU(2) and 
the SU(3) case provided mD was the same jilj. With our choice of an initial anisotropy given by 
to/tiso = 10, equating tq = and using units where tq = 1, the above result corresponds to the 
value m£) = 3.585 employed in Figs. [Mil 

The lower value c ~ 0.5 for the gluon liberation factor corresponds to a smaller Debye mass, 
which turns out to be rather expensive in computer time, because one has then to go to much larger 
values of r to obtain comparable effects and one cannot increase the time steps much without losing 
accuracy. However, in order to see the effect of this lower value of c (which now seems disfavored 
[sHp. it should suffice to simply rescale the r values of Figs. [MSI such that the weak-field Abelian 
regime matches the semi-analytical results presented in Fig. 1 of Ref. [sJ], where c = 0.5 was 
employed. 
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